NAS A/TM— 20 1 0-216961 



Atmospheric Turbulence Modeling for Aero 
Vehicles: Fractional Order Fits 


George Kopasakis 

Glenn Research Center, Cleveland, Ohio 


December 2010 



NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI Program operates under the auspices 
of the Agency Chief Information Officer. It collects, 
organizes, provides for archiving, and disseminates 
NASA’s STI. The NASA STI program provides access 
to the NASA Aeronautics and Space Database and 
its public interface, the NASA Technical Reports 
Server, thus providing one of the largest collections 
of aeronautical and space science STI in the world. 
Results are published in both non-NASA channels 
and by NASA in the NASA STI Report Series, which 
includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or theoretical 
analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 
NASA counterpart of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent of 
graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or 
of specialized interest, e.g., quick release 
reports, working papers, and bibliographies that 
contain minimal annotation. Does not contain 
extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 
http://www.sti. nasa.gov 

• E-mail your question via the Internet to help@ 
sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at 443-757-5803 

• Telephone the NASA STI Help Desk at 
443-757-5802 

• Write to: 

NASA Center for AeroSpace Information (CASI) 
7115 Standard Drive 
Hanover, MD 21076-1320 



NAS A/TM— 20 1 0-216961 



Atmospheric Turbulence Modeling for Aero 
Vehicles: Fractional Order Fits 


George Kopasakis 

Glenn Research Center, Cleveland, Ohio 


National Aeronautics and 
Space Administration 


Glenn Research Center 
Cleveland, Ohio 44135 


December 2010 



Acknowledgments 


The author would like to acknowledge the support of the NASA Supersonics Project for the research conducted in this paper. 


This work was sponsored by the Fundamental Aeronautics Program 
at the NASA Glenn Research Center. 


Level of Review. This material has been technically reviewed by technical management. 


Available from 


NASA Center for Aerospace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 


National Technical Information Service 
5301 Shawnee Road 
Alexandria, VA 22312 


Available electronically at http://gltrs.grc.nasa.gov 



Atmospheric Turbulence Modeling for Aero Vehicles: 

Fractional Order Fits 

George Kopasakis 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


Abstract 

Atmospheric turbulence models are necessary for the design 
of both inlet/engine and flight controls, as well as for studying 
coupling between the propulsion and the vehicle structural 
dynamics for supersonic vehicles. Models based on the 
Kolmogorov spectrum have been previously utilized to model 
atmospheric turbulence. In this paper, a more accurate model 
is developed in its representative fractional order form, typical 
of atmospheric disturbances. This is accomplished by first 
scaling the Kolmogorov spectral to convert them into finite 
energy von Karman forms and then by deriving an explicit 
fractional circuit-filter type analog for this model. This circuit 
model is utilized to develop a generalized formulation in 
frequency domain to approximate the fractional order with the 
products of first order transfer functions, which enables 
accurate time domain simulations. The objective of this work 
is as follows. Given the parameters describing the conditions 
of atmospheric disturbances, and utilizing the derived 
formulations, directly compute the transfer function poles and 
zeros describing these disturbances for acoustic velocity, 
temperature, pressure, and density. Time domain simulations 
of representative atmospheric turbulence can then be 
developed by utilizing these computed transfer functions 
together with the disturbance frequencies of interest. 

Introduction 

This paper addresses the need for a model that simulates 
atmospheric disturbance, in both the time and frequency 
domain, over a wide range of altitudes and variations in 
atmospheric turbulence conditions, that is relatively easy to 
implement and representative of the actual fractional order 
nature of atmospheric turbulence. This is applicable to both 
propulsion system flow field type disturbances as well as 
vehicle gust loads. 

Atmospheric turbulence have been studied for some time, 
especially in the field of Atmospheric Sciences, Nastrom 
(1985), and Fairall (1991), and various models have been 
developed. These models are primarily based on the so called 
Kolmogorov spectrum, originally developed by Tatarski (1961), 
based partly on the studies of turbulence by the Russian 
mathematician Andrei Kolmogorov (Kolmogorov (1941) (a) 
and (b)). The Kolmogorov spectrum has an energy level that 
approaches infinity as frequency approaches zero. This 


characteristic makes it difficult to implement these types of 
models in the time domain. A suitable approximation to the 
Kolmogorov model that is commonly used with a finite energy 
spectrum is the von Karman type model, originally referred to 
as the isotropic-turbulence spectrum, (Houbolt (1964)). 
However, simulating the von Karman type models in the time 
domain is still problematic because of their fractional order. 
The Kolmogorov model has also been extended (Tank (1994)) 
to develop a baseline of atmospheric turbulence for the High 
Speed Civil Transport (HSCT). The Tank model also covers 
atmospheric acoustic wave disturbance modeling utilizing the 
von Karman spectral. Hoblit (1988) introduced the Dryden 
model approximation to the fractional order von Karman 
atmospheric model. But this model is second order compared 
to the 5/3 fractional order of the acoustic velocity atmospheric 
turbulence spectral. Thus, the Dryden model underestimates 
the atmospheric disturbances, increasingly with frequency. 

To alleviate some of these difficulties, the Tank model for the 
von Karman approximation is utilized in this paper to derive an 
explicit electrical circuit analog of atmospheric disturbances. 
This circuit analog will act as a low pass filter. The circuit 
elements are explicitly computed as functions of atmospheric 
parameters, such as eddy dissipation rate and integral length 
scale. However, like the actual atmospheric disturbances, 
Nastrom (1985), the circuit order also turns out to be fractional 
order, which makes it difficult to simulate. Thus, the circuit 
model is used as the basis in this development to derive integer 
order transfer function (TF) approximations to the fractional 
order model. These TF approximations are a product of first 
order poles and zeros, which are determined as a function of the 
parameters describing an atmospheric disturbance. This 
approach alleviates the manual process of hand fitting the 
approximation every time an atmospheric parameter is changed. 

Fairall (1991), Tank ((1994) and (1996)) and others approach 
atmospheric disturbance modeling probabilistically, as an 
exceedance for controls design purposes. That means that the 
probability of a time to failure is computed for controls design. 
This time to failure is associated with a controls design that can 
tolerate a maximum specified disturbance, and computes the 
probability, in terms of flight miles or hours, that an atmospheric 
disturbance will exceed this threshold. In this paper the approach 
will be to compute the worst case disturbance expected, and 
assume that the control system will be designed to handle this 
disturbance. Thus, this paper will not cover exceedance rates , 
which can be computed separately, based on the worst 
disturbance expected and the specifics of the controls design. 
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Kopasakis (2010) extends the work discussed in this paper 
to cover considerations for atmospheric turbulence 
specifications for supersonic propulsion systems. This work 
also provides an example for a supersonic inlet shock position 
controls design in the presence of atmospheric turbulence. 

This paper is organized as follows. First, the Kolmogorov 
form of the atmospheric disturbance spectrum is presented 
followed by the Tank model for the von Karman spectrum 
forms of the acoustic disturbances. This is followed with 
formulations of the equivalent fractional order TF 
approximations of atmospheric turbulence. Finally, time 
domain atmospheric disturbances are discussed for the derived 
TF approximations, followed by concluding remarks. The 
formulations covered in the body of the paper are based on an 
electrical circuit analog developed for different types of 
atmospheric disturbances, which is covered in Appendix B. 
While the detailed derivations of these formulations are 
covered in Appendix C. 


Kolmogorov Form of the Atmospheric 
Disturbance Model 

Tank ((1994) and (1996)) utilizes a Kolmogorov one- 
dimensional locally isotropic atmospheric turbulence 
spectrum, mathematically developed in Tatarski (1961), which 
represents the spectral density of a structured random field of 
atmospheric turbulence as 


S t (k) = a t e 2n k 5/3 ^ 

In this equation the quantity 8 signifies the eddy dissipation 
rate, in units of (energy/(mass x time)) = (m 2 /sec 3 ), and k 
signifies the wavenumber, in units of (rad/m) or (cycles/m). 
The subscript t in the atmospheric turbulence spectral density, 
S t (k ), signifies the type of disturbance. Based on these 
definitions, the units of S t (k) are (m/sec) 2 /(rad/m) (Tatarski 
(1961)), where the units of rad or Hz are not affected by 
raising them to a power. Or S t (k) will have units of 
(m/sec) 2 /(cycle/m) (Tank (1994)), as rad = cycle/2 n, with the 
2tz factor that multiplies the spectral density for this 
conversion absorbed into the constants a t terms. In this 
treatment, more convenient units in terms of Hz will be 
utilized as ((m/sec) 3 /Hz), by substituting cycles with Hz*sec, 
which allows displaying results in terms of the acoustic wave 
velocity in (m/sec) /Hz. 

Tatarski (1961) includes a more detailed treatment of 8 
along with a detailed derivation of Equation (1). In brief, 
s = v * V ^ //^ 

1 ' is called the eddy dissipation rate and it 
represents the energy dissipated as heat per unit mass per unit 
time due to a velocity fluctuation that occurs in an 
atmospheric region of size (length) /, due to a flow instability 
that takes place when a certain critical Reynolds number value 
is exceeded (Tatarski (1961)). This critical Reynolds number 
cannot be determined precisely, but for worst case 


atmospheric turbulence calculations it would be correct to 
assume that this number is exceeded. The quantity v is the 
kinematic viscosity, and v\ is a velocity fluctuation in an 
atmospheric region of size /. 

At lower altitudes, the eddy dissipation rate can vary 
significantly with altitude and severity of turbulence (Fairall 
(1991)) as well as globally (Nastrom (1985)). Tank (1994) 
used a worst value of 8.6e-5 based on data collected at 
altitudes ranging from 25,000 to 40,000 ft (~ 7.5 to 12 km), 
which is about four times the appropriate value of s for North 
Atlantic cruise altitudes. For supersonic cruise altitude, around 
60,000 ft (-18 km), the same value of s can be used, as s is 
relatively constant at altitudes above 20,000 ft (~6 km). 

According to Tatarski (1961), Equation (1) holds for 
turbulence lengths much greater than what would be 
considered a micro structure (10’s of meters) for eddy 
dissipation and much less than the long length scale (that is 
considered to be about 400 km). Some sparse measurements 
performed by Obukhov (1949) in the lower troposphere agree 
with Equation (1), as well as other measurements by Nastrom 
(1985). The constant a t is constant for each type of 
disturbance, given by Tank 1994, to fit observed data: 


a/= 0.15 


a v = 0.2 


ocj 7 — 0.39 

a P = ().m)5(PJT„) 2 


(longitudinal wind velocity gust, 
dimensionless) 

(vertical or horizontal wind velocity 
gust, dimensionless) 

(temperature disturbance, K 2 s 2 m~ 2 ) 
(pressure disturbance, Pa 2 s 2 m 2 ), 
where the subscript “o” denotes total 
atmospheric quantities. 


For a flight vehicle encountering an atmospheric 
disturbance, the disturbance frequency is defined as 


r Ma 7 w 

f = = kMa 

X 


( 2 ) 


Where M is the vehicle Mach number, a is the local speed of 
sound (m/sec), X is the wave length of the atmospheric 
disturbance (m/cycle), and k is l/X (cycles/m). The 
disturbance frequency is based on an aero vehicle traveling at 
a certain altitude and Mach number that encounters an 
atmospheric disturbance with a certain wavelength. There 
would be a corresponding frequency based on Equation (2) 
around which atmospheric disturbances will be exhibited. 
However, a control design would need to be able to 
sufficiently attenuate these disturbances under a wide range of 
flight and atmospheric conditions. Thus, accounting for the 
full operating envelope, a controls design would need to 
attenuate atmospheric disturbances in a wide range of 
frequencies. Normally, lower frequency atmospheric 
disturbances will be easier for the control design to handle. 

The speed of sound is related with temperature as 
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a = yfyRT s 


( 3 ) 


where T s is the static temperature in K. 

According to the International Standard Atmosphere (ISA) 
or (Anderson (2000)), the temperature decreases at a rate of 
approximately 6.5 K per km, from sea level up to 36,000 ft 
(-11 km). From 36,000 to 65,000 ft (-11 to 20 km) which 
would be the approximate cruise altitude of a supersonic 
transport, the temperature remains constant at about 216 K. 
Thus, temperature (K) as a function of altitude, h (km), can be 
expressed as 


T(h) = T oh -6.5(h-h oh ) (4) 

where r oh is the temperature at an arbitrary reference height 
/z oh , which is accurate up to 36,000 ft and then stays constant 
above that. Based on that, at supersonic cruise conditions the 
speed of sound would be approximately 295 m/sec. 

Figure 1 shows the zonal (east- west) winds, meridional 
(north- south) wind, and potential temperature power spectral 
as reported in Nastrom (1985). The units of the wind power 
spectral are in (m 3 /sec 2 )/rad and temperature is in (K 2 m/rad). 
As shown in this figure, at the long length scale (i.e., beyond 
400 km) the power of turbulence decreases as the -3 power 
with decreasing wavelength (increasing frequency), thereafter 
decreasing as the -5/3 power which agrees with Equation (1). 
Therefore, Equation (1) shows a good representation of the 
atmospheric spectral density with a wavelength smaller than 
400 km. 



Figure 1. — Wind and Potential Temperature Spectra as 
reported by Nastrom (1985). Note: for clarity, the meridional 
wind and potential temperature spectra have been shifted 
one and two decades to the right, respectively. 


The Tank Model of the Atmospheric 
Disturbance 


As an alternative to the Kolmogorov spectrum, Tank 
(Soreide and Tank (1996) and (1997)) scaled the von Karman 
spectrum to fit the Kolmogorov model in the limit (i.e., for 
large k ), which compares well with other data. For 
longitudinal disturbances this spectrum is 


Si,VK:{k)= 2.7s 2/3 i 5/3 


2 

l + (l .339(271 )i&) 2 


( 5 ) 


For transverse disturbances this spectrum is 


Sv,VK{k) 


\ + -(l.339(2n)Lkf 

2.7s 2/: V 5/3 2 

l + (l.339(27t)Zi:) 2 ] 11/6 


( 6 ) 


The main difference of this von Karman form (compared to 
Kolmogorov model) is that the disturbance spectral levels off at 
low frequencies. Another difference is that in the von Karman 
form of the Tank model, the integral length scale, Z, is explicitly 
employed, Equations (5) to (6). The integral length scale is 
related to the outer scale length, L 0 (the length of the 
atmospheric turbulence patch) by L=LJ\A.l for the longitudinal 
disturbance, and L=2LJ\9 for the transverse disturbance 
(Soreide and Tank (1996)). Figure 2 shows a plot of the 
Kolmogorov acoustic wave velocity spectral in (m/sec/Hz) 
compared to the corresponding von Karman spectral for the 
longitudinal and transverse acoustic wave velocity disturbances. 
Different values of L will produce the same curves 



Figure 2. — Acoustic wave velocity spectral comparisons for the 
Kolmogorov and von Karman spectral. 
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Frequency, Hz 

Figure 3. — Longitudinal acoustic wave velocity spectral 
comparisons between the Kolmogorov and the von Karman 
forms of four different integral scale lengths. 



Frequency, Hz 

Figure 4. — Longitudinal acoustic wave velocity spectral 
comparisons for different integral scale lengths. 


with either a higher or lower, low-frequency asymptote as 
shown in Figure 3 for the longitudinal disturbance. For 
comparison, a value of Z=762 m pertains to an atmospheric 
turbulence patch of approximately 1 1 km for the longitudinal 
disturbance, while a value of Z=10 km pertains to a turbulence 
patch of approximately 147 km long. Tank used a value of 
L=162 m, which he called as standard in the airplane industry. 
Since a typical control design can sufficiently attenuate 
disturbance frequencies well beyond 1 Hz, based on Figure 3 
differences in the value of L and thereby, differences in the 
lower frequency asymptote, will have negligible effect on the 
control design. For completeness, the Kolmogorov spectrum 


for the longitudinal acoustic wave velocity is shown in 
Figure 4 for different values of eddy dissipation rates. As can 
be deduced from Figure 4, orders of magnitude difference in 
eddy dissipation rate doesn’t produce a corresponding 
appreciable change in the Kolmogorov spectral density. For 
instance, approximately four orders difference in s only makes 
approximately a factor of 8 difference in velocity as 


W) = 


/ x2/9 

£2 


V c l J 


S\(f) , ((m/sec)/Hz) 


( 7 ) 


Where S\ is a known spectral density with an eddy dissipation 
rate c i , and S 2 is the calculated spectra density for a different 
eddy dissipation rate s 2 . 


Fractional Order Fit of Atmospheric 
Turbulence Model 

Atmospheric turbulence, as shown in Figure 1, and as 
described by Kolmogorov and the Tank models in 
Equation (1) and Equations (5) and (6), is fractional order. In 
Appendix B a circuit analog of the Tank model von Karman 
form is utilized, which serves as the basis for deriving integer 
order pole-zero product TF approximations to the fractional 
order atmospheric disturbances. The reason for the need to 
derive approximations to the fractional order equations is 
because of the difficulty of explicitly or numerically solving 
fractional order differential equations. The reason is fractional 
order derivatives, unlike integer order derivatives, do not obey 
the locality law (i.e., the limit theorem), as further explained in 
Appendix B. 

The idea behind the integer order TF approximation for a 
fractional order TF is as follows (see Fig. C.l): Starting at a 
frequency near the beginning of the equivalent 3 dB point of 
the fractional order TF, an integer order TF approximation can 
be developed symmetrically centered about the fractional 
order TF, like a descending staircase shape, by interleafmg 
integer order poles and zeros. As the number of steps of this 
staircase TF approximation increase, the steps become shorter, 
eventually collapsing to the straight line of the fractional order 
TF as the number of poles and zeros of this approximation is 
increased to infinity. For this approach to work, the frequency 
of the poles and zeros need to be related to the atmospheric 
disturbance parameters and also be derived in a way such that 
the staircase TF approximation is symmetrically centered 
about the fractional order TF. 

Based on the fractional order circuit analog of Equations (5) 
and (6), see Appendix B for detailed derivations, the values 
for the equivalent capacitance and resistance of atmospheric 
turbulences are determined as 


(^s 2/3 ) 1/x (27iMa) 


( 8 ) 
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R, = 1.3390(2n)(a,e 2/3 f x L 


( 9 ) 


and the corresponding natural frequency is computed as 


m z = m p - 1 (16) 

For a desirable pole-zero density pair, p pz , per decade to be 
used to approximate the fractional order disturbance 


0) 


n 


Kg m 

m 


(10) 


1 

2ppz 


(17) 


where the adjustment factor K (on is 1, but represented in this 
form in case any adjustments need to be made to the 
atmospheric disturbance natural frequency. 

Utilizing these circuit analog parameters, a time domain 
disturbance can be derived (see Appendix C for detailed 
derivations) based on an integer order TF approximation for 
atmospheric turbulence, given the atmospheric disturbance 
parameters s, L , q , t , and r for the units conversion factor (with 
q=xr ; x=5/3 and 7=1/3 for acoustic disturbances and Vi for 
temperature and pressure) as 

m z 

P[(.v/o) r , +])) 

W t , 0 =K tJlt - L W t (11) 

n(v°v+i) 

1 

where W t is series of input sinusoids with unit amplitude 
frequency components distributed over the frequency range of 
interest, which is discussed more in Kopasakis (2010) for 
propulsion system type disturbances. The frequencies of the 
poles can be computed as 

i — 1 

Kapi&Hpi ( ®Hpi /(fipi-j + l) 

( 0 pi = i = 2,3,...,m p (12) 

1 0 r|(2i- 1 ) 9 ]q (o:> IIp i joi z j_j + 1 )- 1 

7=1 


and the terms (o Hpi and <x> Hzi in Equations (11) and (12) can be 
computed as 

®flp/=® B (l0 119(2M) -lf 9 5 i = 2,...,m p (18) 

m H Z i =oj„(l0 2 'W' ~lf ,q , i = \,2,...,m z (19) 

The utility of (o Hpi and co Hzi are to maintain symmetry, so that 
the staircase pole-zero approximation is symmetrically located 
on top of the fractional order TF asymptote. The proportional 
gains K (opi and K (0zi in Equations (12) and (14) are reserved for 
final adjustments that may be needed to these frequencies for 
more closely approximating the fractional order TFs 
representing atmospheric disturbance. 

Longitudinal and Transverse Acoustic Wave 
Turbulence 

The longitudinal and transverse acoustic wave atmospheric 
disturbances are in the form of pure wind gusts in the axial 
direction of vehicle motion and in the vertical direction of 
motion respectively. Based on the fractional order fits 
determined in Appendix C, using Equations (8) to (19), with 
q=xr= 5/9 (r=l/3), then K t> f lt in Equation (11), based on 
inspection of Equations (5) and (6), is 

Kifit =(5.4 £ 2 / 3 jL 5/3) 1/ 3 (20) 


where the first pole is computed as follows 


for the longitudinal disturbance, and 


(tipi 



(13) 


^v,fi t =(2.7s 2/3 I 5/3 ) 1/3 (21) 


and the frequencies of the zeros can be computed as 


K($zi^Hzi Y | i^Hzi /®. zi-j l) 

co zi = i=l,2,...,m z (14) 

10" 2 w ^[(co^/co^ +l)-l 

7=1 


for the transverse disturbance. For n = 3 (i.e., the TF fit over the 
span of 3 decades in frequency), the proportionality factor 
adjustments to improve these fits have been found to be 
(Appendix C), 

*/,<*> = [^co n ? ^ pi , K($zi ] (22) 

= [2.4; 1 1 1/2.4 1/1.5; 1 1 l] 


For n number of frequency decades desired to be estimated, the 
number of poles and zeros in the TF approximation can be 
determined as 


m p = (n- 1 )/q 


(15) 


for the longitudinal acoustic wave, and 


K 


v,co " 


0 pi , J 

[4.27; 1 1 1/2.4 1/1.5; 1 1 l] 


(23) 
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for the transverse. These TF fits are compared against the 
Tank von Karman spectral of Equations (5) and (6), with the 
units conversion exponent r, as shown in Figures 5 and 6 for 
s=8.6e-5 (m 2 /sec 3 ) and L=162 m. As shown in these figures, 
the fits do a good job in approximating these disturbances. 
More accuracy can be achieved around the knee of these 
spectral by increasing the density of pole-zero pairs per 
frequency decade in this region. A certain inaccuracy at the 
lower frequencies would be acceptable for a feedback control 
design, since a typical control design should have no problem 
attenuating disturbances at this low frequency range. 



Figure 5. — Longitudinal von Karman spectral, final adjusted TF 
fit (e = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Figure 6. — Transverse von Karman spectral, final adjusted TF 
fit (e = 8.6e-5 m 2 /sec 3 , L = 762 m). 


Temperature Turbulence 

Atmospheric temperature turbulence causes both temperature 
as well as acoustic velocity disturbances due to the change in 
the speed of sound. For the most part, temperature disturbances 
result in vertical displacement of air parcels (so called gravity 
waves). Therefore, acoustic velocity disturbances caused by 
temperature will generate vertical wind gusts that add with any 
transverse acoustic velocity gusts. However, for the propulsion 
system, according to Ashun (2004), the wing forebody will turn 
a vertical gust into a longitudinal gust, by multiplying the 

vertical gust by the conversion factor ( M ^ - 1)/ -1 . This 
factor amounts to 0.63 for M ^ = 2.35, and can be ignored for 
worst case purposes. However, this conversion factor 
decreases with speed and can be taken into account, especially 
at lower supersonic speeds, for longitudinal acoustic gusts due 
to temperature turbulence. In Appendix B, Section B.2, first the 
Kolmogorov spectrum of the longitudinal acoustic wave and 
that of the temperature, based on Equation (1) are plotted 
(Fig. B.6). This results in parallel spectral lines with frequency, 
similar to the Kolmogorov spectra shown in Figure 2. Then the 
von Karman temperature spectral is constructed by scaling the 
horizontal von Karman type acoustic wave spectral (i.e., scaling 
it by the difference in magnitudes of these Kolmogorov spectra 
as shown in Fig. B.6). 

Temperature turbulence spectral density, as can be seen in 
Figure 1, also follows the 5/3 law, but its units are in terms of 
Kelvin squared. Thus, in order to convert the units to Kelvin, 
the exponent r becomes 14, which makes the fractional order 
q=5/6. Based on this, the TF fit performed for temperature is 
done the same way as with acoustic wave in the previous 
section (i.e., utilizing Equations (8) to (19), with the additional 
relations. 

^,f,t(temp)-Vl4.0s 2/3 L5/3 (24) 

Kt,g> = Kn > K (1 )pi , K azi ] (25) 

= [l.5; 1 1 1/1.1 1/1.2; 1 1 l] 

Figure 7 shows this TF fit and the spectral of Equation (5) 
(i.e., Eq. (5) substituting the scaled Eq. (24) for the numerator 
and by also applying the units conversion factor in the 
denominator, for s=8.6e-5 anZ=762 m. 

The proportionality factors, K h for all the fits can also be 
considered part of the formulations. Their accuracy at higher 
frequencies starts deviating somewhat with higher values of L 
(integral length scales of few thousand meters and above), 
values that are considered to be outside the typical range. As 
discussed before, Z=762 m is considered standard in airplane 
industry. 

Temperature turbulence also causes acoustic wave gusts due 
to change of the local speed of sound. By perturbing the 
relation of the speed of sound and temperature in Equation (3), 
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Figure 7. — Temperature von Karman spectral and 
its TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


and substituting the resulting expression for A a into the Mach 
number equation of M=v/a , the relation between the change in 
temperature and acoustic velocity can be obtained as 


Av = 


My R ArT7 
— *—AT 
2u 0 


(26) 


This relation can be applied to Equation (24) to compute 
^r,fit(acoustic) for the acoustic wave velocity disturbance due to 
an atmospheric temperature fluctuation as 

^7\fit( acoustic) = ^Vl4.08 2/3 I 5/3 (27) 

Z(2 0 

Equation (27) is applicable to transverse or vertical wind gust 
disturbances generated by temperature turbulence, and 
accurate in the worst case sense for longitudinal wind gusts, 
especially at higher supersonic Mach numbers. Generally, 
longitudinal wind gusts due to temperature variations can be 
expressed as follows by utilizing the conversion factor 

(M 0 o -1 )/Vm 2 -1 

M x {m x -l)y 7 ? j „ 2/3 7-5/3 po\ 

A 7\fit(longitudinal) - VI 4.1)8 L (28) 

2 a 0 ^JMl -1 

This TF fit for the acoustic velocity disturbance at supersonic 
cruise altitude, at Mach 2.35, for L=762 m and two different 
values of 8 is shown in Figure 8. The TF fit for e=8.6><10” 5 
(m 2 /sec 3 ) and for two different values of integral scale lengths 
is shown in Figure 9. As can be seen from Figure 8 or Figure 9 
compared to Figure 5, the acoustic velocity gusts produced by 
fluctuations in temperature are much higher in amplitude than 
those produced by pure acoustic velocity gusts. 



Figure 8. — Von Karman acoustic wave velocity spectral due to 
temperature gust and its TF fits for different eddy dissipation 
rates (L = 762 m). 



Frequency, Hz 

Figure 9. — Von Karman acoustic wave velocity spectral due to 
temperature gust and its TF fits for different integral scale 
lengths (e = 8.6e-5 m 2 /sec 3 ). 


As shown in Appendix B, Figure B.8, for the combined 
acoustic wave, the acoustic wave velocity due to temperature 
gust dominates in the low frequency range and at higher 
frequencies both acoustic waves affect the total. At worst case, 
it can be assumed that both affects of the acoustic wave 
velocity can combine. It is also assumed that both the 
longitudinal and transverse acoustic wave velocities will not 
combine to produce worst case affects. For engines under the 
wing, the transverse wave is assumed to be converted to a 
longitudinal wave via the wing forebody (Ahsun (2004)), see 
Appendix B. 
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Pressure Turbulence 


The longitudinal velocity wave spectral and the pressure 
spectral of Equation (1) was utilized the same way as 
described in the previous section for the temperature to scale 
the von Karman type acoustic wave spectral to come up with 
the pressure spectral of atmospheric disturbances, see 
Section B.3. The units conversion exponent for pressure is the 
same as that for temperature (i.e., r= 1/2), which makes the 
fractional order the same (i.e., q= 5/6). As such, K pflt and K p & 
were computed as follows 

K P , flt = Vll.6s 2/3 L 5/3 (29) 

Kp,® — ho m > K a pi , K az i ] (30) 

= [l .5; 1 1 1/1.1 1/1.2; 1 1 l] 

Figure 10 shows this TF fit for e=8.6e-5 (m 2 /sec 3 ) and 
L=162 m compared to the scaled von Karman type spectral — 
i.e., Equation (29) substituted for the numerator and 
proportional factor of Equation (5) and the denominator of 
Equation (5), raised to the power r. 


Density Turbulence 

Density disturbances are contained within temperature and 
pressure fluctuations and its inclusion in a simulation, along 
with pressure and temperature will produce unnecessary 
additional disturbances. However, if need be, the relation of 
density disturbance with temperature and pressure can be 
derived as (by perturbing the state equation) 


Ap = P ° +AP -A. 

R(T 0 +AT ) RT a 


(31) 


A plot of density disturbance utilizing Equation (31) is 
shown in Figure 11. 


Simplified Atmospheric Turbulence Model for 
Controls Design 

The attempt in this section will be to further simplify the 
fractional order TF fit developed in the previous sections, by 
identifying certain characteristics of this model that are 
important for controls design to handle atmospheric 
turbulence. 

Based on observations, utilizing Figure 8, differences in 
eddy dissipation rates, 8, produce spectral densities with the 
same frequencies, offset by a fixed magnitude. This offset is 
accounted by differences in TF gain. Calculating these 
frequencies using the formulations developed so far also 




Frequency, Hz 

Figure 1 1 . — Density von Karman spectral and its TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


supports this observation. Differences in the fractional order of 
the disturbances, such as longitudinal acoustic wave velocity 
vs. temperature, obviously, will produce different spectral 
frequencies. Unlike differences in eddy dissipation rates that 
produce the same frequencies with different TF magnitudes, 
inspection of Figure 9 shows that spectral with different 
integral length scales, L , will be represented by different 
frequencies. However, different integral length scales only 
affect the low frequency spectral, where a typical feedback 
controls design should have no difficulty attenuating 
disturbances at such a low frequency range. 
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Figure 12 shows spectral density TF fits, first, with a fixed 
integral length scale and two different values of 8. This shows 
that the differences in s produce a fixed offset in magnitude, 
but with the same pole-zero frequencies. Then Figure 12 also 
shows comparison between two spectral densities with the 
same value of 8 and two different values of L. This later 
comparison, with different values of Z, is represented with 
different pole-zero frequencies. However, their high frequency 
asymptotes match, which means that for the purpose of 
controls design, either of these TF fits can be utilized, despite 
their differences in pole-zero frequencies. 

Therefore, the frequency domain atmospheric turbulence 
model can be simplified by utilizing the model developed in 
the previous section to calculate the fixed frequencies due to a 
fixed integral length scale and variable eddy dissipation rates. 
Based on that, the simplified model developed for atmospheric 
turbulence, which ignores differences in the pole-zero 
frequencies due to different integral length scales can be 
represented as 


Where G LA , G VA , G T , G TLA , G TVA and G P are the simplified 
atmospheric disturbance TF for longitudinal and transverse 
acoustic wave velocity, temperature, longitudinal and vertical 
acoustic velocities due to temperature, and pressure 

respectively. The factor - 1 )/^M£ - 1 in Equation (35) 
can be set to one for worst case turbulence conditions, 
especially for relatively high Mach numbers, in which case it 
makes Equation (35) and Equation (36) the same. 

The TF fit described in the previous section was carried up 
to 200 Hz, which should make it accurate even beyond this 
frequency range. However, the model simplification 
represented by Equations (30) to (34) introduces additional 
error at high frequencies and its validity should be limited to 
about 200 Hz. The frequencies of this model are based on an 



Figure 12. — von Karman longitudinal acoustic wave velocity 
spectral and its TF Fits due to different values of eddy 
dissipation rates and Integral Length Scales. 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 

integral length scale of L = 762 m. This means that this 
simplified model is not unique, and other equivalent models 
can be calculated using different values of L with the 
equations developed in the previous sections. Thus, for a 
controls design, which can adequately attenuate 1 Hz 
disturbances, this model is accurate for values of L as low as 
100 m, which is well below the expected range of L. A value 
of L=162 m is considered standard in the airplane industry as 
mentioned before. 

Time Domain Disturbances 

The TF fits in the previous section can be thought as filters, 
which take input disturbances at desirable frequencies and 


q (s) 10£ 2 /9 (W9.2 + 1)(W55.0 + 1)(W335.5 + 1) 

U[ ’ (5 / 1 .46 + 1)(5 / 30. 1 + l)(s / 85.7 + 1)(5 / 1 593. 1 + 1) 

Gy A (s) = 5682/9 (W9.2 + l)(W55.0 + l)(,/335.5 + l) 

(5 / 1 .46 + l)(s / 30. 1 + l)(s / 85.7 + l)(s / 1593. 1 + 1) 

q (s) = 943£ 2/6 ( 5 / 33.0 + !)($/ 45.6 + !)($/ 602.4 + 1) 

(s/1.1 + 1)(5/25.1 + 1)(s/109.8 + 1)(5/816.3 + 1) 

G (j) ms 2/6 (M ~ 1} (j/33.0 + 1)(j/45.6 + 1)(j/602.4 + 1) 

TLAK ’ Vm 2 - 1 a Q ( 5 /I.I + 1)(5 / 25. 1 + l)(s / 109.8 + l)(s / 8 1 6.3 + 1) 

G (~) 172c Cs733.0 + l)(s745.6 + l)(s7 602.4 + 1) 

TVAK ’ a 0 ( 5 /I.I + 1)(5 / 25. 1 + 1)(5 / 1 09.8 + 1)(5 / 8 16.3 + 1) 

G ( 5 ) g5?c 2/fi (W 33.0 + l)(5/45.6 + l)(5/602.4 + l) 

M J ( 5 /I.I + 1)(5 / 25. 1 + 1)(5 / 1 09.8 + 1)(5 / 8 1 6.3 + 1) 
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convert them to representative free stream atmospheric 
disturbances. The TF fits contain the representative 
magnitudes of these disturbances. Thus, the discrete time 
domain frequencies, W h only need to have unity RMS values. 
For instance, a time domain input can be constructed as a sum 
of unit amplitude sinusoids, starting with a low frequency, 
somewhere at the low frequency asymptote of the spectral, 
and continuing up with sinusoids that are dispersed along the 
high frequency asymptote up to frequencies that covers the 
controller bandwidth. Also, the control design can be checked 
with single frequency sinusoids, sequentially simulated. Of 
course, the total energy or the maximum amplitude of the 
disturbance will need to be reasonable. For instance, it may be 
considered that total wind gust will not exceed an amplitude of 
180 mph, or some other reasonable number that depends on 
altitude. For more information on propulsion type disturbance 
considerations, see Kopasakis (2010). 

For most of the spectral density simulations covered so far, 
the worst case eddy dissipation rate for low to moderate 
turbulence was used for altitudes of about 36,000 ft (-11 km) 
and above. However, the highest atmospheric turbulence ever 
recorded had an eddy dissipation rate of 1.7e-3 (m 2 /sec 3 ) 
(Tank (1994)). Theoretically, it is possible that even more 
severe turbulence can be encountered according to McMinn 
(1997). McMinn provides tables for altitudes, turbulence 
severity, and associated eddy dissipation rates. Kopasakis 
(2010) covers more in this area of simulating atmospheric 
disturbances, including considerations for atmospheric 
turbulence specifications. Fairall (1991) also covers 
implementation of time domain sinusoidal input disturbances 
by utilizing Fourier sinusoids. 

As mentioned before, previous works cover atmospheric 
turbulence probabilistically, as an exceedance for controls 
design purposes. In terms of probabilistic atmospheric 
turbulence models, there are no comparisons drawn here 
between those models and the non-probabilistic approach 
described in this paper, even though exceedance can also be 
incorporated with this approach. The modeling approach 
developed in this paper is geared towards implementing 
controls design by considering worst case atmospheric 
turbulence. This also allows implementation of control designs 
to cover worst case turbulence encounters, without penalizing 
performance under normal atmospheric conditions, by 
incorporating control scheduling in the design. 

Figure 13 shows a time domain longitudinal acoustic wave 
disturbance utilizing the acoustic wave TF fit discussed in the 
previous section. The input to this TF is a combination of 
unity amplitude sinusoids distributed along the frequency 
range shown in Figure 5. For a single unity amplitude sinusoid 
at low frequency, the amplitude of the acoustic disturbance 
would be approximately 9 m/sec according to Figure 5. For 
this summation of unity amplitude input sinusoids, their 
amplitudes combine to give approximately 15 m/sec 



Time, sec 

Figure 13. — Longitudinal acoustic wave velocity disturbance 
created using the TF fit of the longitudinal acoustic wave 
with unity input sinusoids. 


maximum disturbance. The same approach can be used to 
construct other atmospheric disturbances for time domain 
simulations. 

In Kopasakis (2010), considerations for atmospheric 
turbulence specifications are provided, with more discussion 
on the frequencies of these turbulences, also involving the 
vehicle Mach number. In this reference a control design 
example is also provided for the supersonic inlet shock control 
problem, in the presence of atmospheric turbulence. 

Conclusion 

In this paper derivations were carried out to approximate the 
fractional order atmospheric turbulence with integer order 
pole-zero transfer function fits for more accurate time domain 
simulations. An existing von Karman type model form of the 
Kolmogorov spectral is utilized initially for the acoustic 
velocity gusts disturbances. This model is scaled in this paper 
in order to also develop the von Karman model forms for other 
type of atmospheric disturbances like temperature, pressure, 
and density. Because of the fractional order of these models, a 
circuit equivalent is developed that is used as the basis to 
derive the integer order pole-zero approximations. Utilizing 
the formulations presented in this paper, the transfer function 
poles and zeros of the approximations to the atmospheric 
disturbances can be directly computed for different parameters 
describing atmospheric turbulence. The model is derived to 
duplicate the fractional order form of atmospheric turbulence 
and is more accurate than previous models developed. These 
new models could be used to design controls for aerospace 
vehicle propulsion or flight control systems. 
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a speed of sound, (m/sec) 

a Q ambient speed of sound, (m/sec) 

C t equivalent electric circuit capacitance for t-type of 
disturbance, (farads) 

/ frequency, (Hz) 

f n natural frequency, (Hz) 

f c frequency of computed correction factor, (Hz) 

h altitude, (km) 

H vector of intersection points to match estimated and 
fractional order TFs 

k wavenumber, (rad/m) or (cycles/m) 

K proportional gains 

K t von Karman horizontal asymptote for t-type of 
disturbance 

K an adjustment factor to disturbance natural frequency 
/ length, (m) 

L integral length scale, (m) 

L 0 outer scale length, (m) 

m p number of poles in the TF approximation 

m z number of zeros in the TF approximation 

M mach number 

n number of frequency decades desired to estimate the 

fractional order TF 

P pressure, (Pa) 

P Q standard atmospheric pressure, (Pa) 

q fractional order of equivalent electrical circuit 

R universal gas constant, 287 (N*m)/(kg*K) 

r units conversion exponent of atmospheric disturbance 

R t equivalent electric circuit resistance for t-type of 
disturbance, (ohms) 

S atmospheric turbulence spectral density 

Si- or-v atmospheric velocity turbulence spectral density, 

(m 3 /sec 2 )/(rad) or (m/sec) 3 /(Hz) also converted to 
(m/sec)/(Hz) by taking 1/3 root 

s Laplace operator 

T temperature, (K) 

T 0 standard atmospheric temperature, (K) 

T oh temperature at an arbitrary reference height, (K) 
v flow velocity, (m/sec) 

v\ velocity fluctuation in a region of size / of the basic 

laminar flow, (m/sec) 

W unity disturbance time domain signal 

W t disturbance time domain signal for t-type of 

disturbance 


fractional order of atmospheric disturbance spectral 

A.l Greek 

a constant associated with atmospheric turbulence 

spectral density 

a / longitudinal, (unit less), 

a v transverse or vertical, (unit less) 

a T temperature, (°K 2 sec 2 m" 2 ) 

clp pressure, (Pa 2 sec 2 m" 2 ) 

y ratio of specific heats, (y = 1.4 1) 

8 eddy dissipation rate, (m 2 /sec 3 ) 

r| the ratio of each decade interval where a pole or a 

zero will be used to estimate the fractional order TF 

L wave length of atmospheric disturbance (m/cycle) 

v kinematic viscosity of a laminar flow, (m 2 /sec) 

p weight density, (kg/m 3 ) 

p pz density of pole-zero pairs per decade for estimated 

TF 

co„ natural frequency, (rad/sec) 

(dpi frequency of TF pole i of integer order approx, 
(rad/sec) 

co zz frequency of TF zero i of integer order approx, 
(rad/sec) 

A.2 Subscripts 

a constant associated with circuit capacitance 

correctional factor 

A variable associated with an adjustment 

C constant associated with circuit capacitance 

correctional factor 

H associated with symmetry frequencies for the approx, 
fit variable associated with TF fit or approximation 

/ variable associated with longitudinal atmospheric 

disturbance 

LA variable associated with longitudinal acoustic 

disturbance 

VA variable associated with vertical or transverse 

acoustic disturbance 

P variable associated with pressure disturbance 

t variable associated with type of disturbance 

T variable associated with temperature disturbance 

TLA variable associated with longitudinal acoustic 

disturbance due to temperature 


Appendix A. — Nomenclature 

X 
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TVA 

v 

VK 

VKA 


variable associated with vertical or transverse 

acoustic disturbance due to temperature 

variable associated with vertical or transverse 

atmospheric disturbance 

variable associated with the von Karman spectral 
density 


variable associated with the von Karman spectral 
density circuit approximation 


o 

P 

s 

z 

(Dpi 

(Dzi 

00 


variable associated with an output 
variable associated with TF poles 
variable associated with a static quantity 
variable associated with TF zeros 
variable associated with pole frequencies 
variable associated with zero frequencies 
variable associated with freestream conditions 
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Appendix B. — Circuit Analog 

In this appendix the von Karman approximations to the 
different atmospheric disturbances based on the Kolmogorov 
spectrum given by (Soreide and Tank (a) and (b)) are 
approximated by utilizing an electrical circuit analog. The 
purpose of the circuit analog is to serve as basis for deriving 
integer order pole-zero product TF approximations to the 
fractional order atmospheric disturbances. The reason for the 
need to derive approximations to the fractional order equations 
is because of the difficulty of numerically solving fractional 
order differential equations. For instance, unlike integer order, 
fractional order derivatives obey the law of non-locality. This 
property of fractional order derivatives makes the state 
transition matrix a convolution integral, with each integration 
step necessitating the computation of the convolution of the 
state all the way back to time zero where the state had an 
initial value of zero. 


B.l Longitudinal and Transverse Disturbances 

The Kolmogorov atmospheric turbulence spectral model as 
described in the main body, Equation (1), of this paper is 

S t (k) = a t e 2/3 k ~ 5/ 3 ( B j) 

And the von Karman longitudinal and transverse models 
respectively, reproduced from Equations (5) and (6) from the 
main body of this report are: 


S, yK {k) = 2.7s 2li L 5/i 
Sv,VK:{k) = 2.7s 2/3 Z, 5 /3 


2 

l + (l.339(2?r)«) 2 5/6 

1 + |(1.339(2ti)Z,A:) 2 
1 + (l . 339(271 A:) 2 U/6 


(B-2) 


(B-3) 


The shape of the von Karman spectrum approximations of 
Equations (B.2) and (B.3) is that of a parallel network 
consisting of a resistive element corresponding to the horizontal 
von Karman asymptote and capacitive element corresponding to 
the high frequency Kolmogorov asymptote as 

S t ,VKA {k) = ; (m/sec)/Hz (B.4) 

{K t ) r +S[{k) 

Where r represents the units conversion exponent from 
(m/sec) 3 /Hz to (m/sec)/Hz (i.e., r=l/3 in this case — for the 
longitudinal and transverse acoustic wave disturbances). The 
variable S t (k) signifies the Kolmogorov spectrum of 
Equation (B.l) or the high-frequency asymptote of the von 
Karman spectrum, and K t signifies the von Karman low- 
frequency asymptote or TF gain for the type of disturbance of 
Equations (B.2) and (B.3) as 


of the von Karman Model Form 


K, =5.4e 2/3 L 5/3 
K v =2.7e 2li L 5n 


(B.5) 

(B-6) 


Equation (B.4), for a flow, such as that of the acoustic velocity 
spectral, is equivalent to a circuit as shown in Figure B.l. The K t 
element is a gain for the type of disturbance and is a function of 
the parameters describing the atmospheric disturbance. The 
parameter W t represents an input flow and W Uo represents the 
resulting output flow, where the exponent q represents the 
fractional order of these elements. The circuit is an analog of 
this representation, acting as a disturbance source connected to a 
filter, meant to be connected to a sink or a system simulation in 
order to establish an output flow. Applying circuit theory, the 
TF of this circuit can be computed as 


W t , 


K[ 


{RtCtsf+l 


W t 


(B.7) 


Where W t in the circuit stands for a unity flow signal 
representing the atmospheric disturbance, containing the 
disturbance frequencies of interest. More on that will be 
discussed later. The output signal, W t>0 , will be the time 
domain signal representing the von Karman spectral 
approximation. The subscript t stands for the type of 
disturbance as before, the exponent q represents the fractional 
order of the circuit, and the factor K t is given by Equation (20) 
or (21). Where l/(C , t s , ) q and R t q are the capacitive and resistive 
impedances of the circuit, with s representing the Laplace 
operator. 

The values of resistance R h capacitance C t , and the 
fractional exponent q remain to be determined. The natural 
frequency of this circuit, co„, can be directly obtained by 
inspection of the TF of Equation (B.7). This frequency will 
coincide with the frequency at which the low-frequency 
asymptote (Figs. 2 and 3) intersects the high-frequency 
asymptote associated with the Kolmogorov spectral or the 
capacitive impedance of the circuit as 



R, q 

vwv — 0 


» w u 



1/(C,s)F 


Figure B.l . — Equivalent circuit TF of atmospheric 
disturbance. 
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CO. 


(B.14) 


For now K^ n is used as a placeholder in case any adjustments 
are needed to this natural frequency and its value here is set to 1 . 

In Equation (B.l), k stands for the wavenumber and it is 
related to the frequency of the disturbance as 

f = = (B.9) 

X 

For the circuit analogy, the natural frequency of the circuit 
should be set to the same as that for the von Karman spectral 
of Equations (B.2 and B.3). By inspection, the natural 
frequencies of both the von Karman longitudinal and 
transverse disturbances are the same and can be calculated by 
substituting Equation (B.9) for / into Equation (B.2) or (B.3) 
to obtain the relationship 

— = l.339(2n)L^-. (B.10) 

co v ’ Ma 

Then by substituting the magnitude s for 2nf on the right hand 
side of Equation (B.10) and also setting co„ = 2nf n on the left 
hand side, we can solve for the natural frequency as 


Additionally, the TF of Equation (B.7) needs to represent the 
von Karman approximation of the Kolomogorov model (i.e., 
the Kolmogorov asymptote). However, the Kolomogorov 
model, as a straight line decreasing with frequency in a log 
scale, can be represented in Figure B.l with a capacitive 
impedance. Thus, Equation (B.l) must be this capacitive 
impedance. This capacitive impedance is raised to the r=l/3 
power to account for the units change from (m 3 /(sec 2 cycle)) to 
(m/sec)/Hz. Therefore, taking Equation (B.l) and raising it to 
the power r, replacing k with f/{Ma) using Equation (B.9), then 
multiplying both numerator and denominator by (27i) (5/3)r , and 
by setting the resulting equation to l/(C^) g , the following 
expression can be obtained. 


By inspection of Equation (B.12), since |s| = 2nf, then 
q= (5/3 )r or in general q=xr (i.e., with x=5/3 in this case), then 
the capacitance, C h can be computed as 

C t —~ (B.13) 

(a t z 2/3 J x (2nMa) 

Plugging Equations (B.ll) and (B.13) into (B.8), with 
co = 2iif R t can be solved as 


R, =1.339(2tc fa t £ 2/3 J /x L 

Neither the capacitance C t nor the resistance R t depends on the 
exponent r, utilized for units conversion, which could be 
indicative of the fundamental nature of these parameters. 

With this analysis, the circuit elements for the longitudinal 
and transverse von Karman model approximations, pertaining 
to Equation (B.7) have been computed. Figures B.2 and B.3 
show the von Karman longitudinal and transverse spectral of 
Equations (B.2) and (B.3) respectively and their circuit 
approximation of Equation (B.7). These approximations can 
be improved by increasing the magnitude of the circuit 
approximation (K t in Equation (B.7) for the longitudinal wave 
in order to reduce the high frequency asymptote error at the 
lower frequency. 

To do this, an offset dB level is chosen (AdB) so that the 
resultant horizontal asymptote of the circuit approximation is 
not too high above the corresponding horizontal asymptote of 
the von Karman spectral. This can be done by multiplying 
Equation (B.7) by a proportional gain K a , corresponding to 
this AdB level, as K a =10 AdB /2 °. However, doing that will also 
raise the high frequency asymptote of the circuit 
approximation by the same dB level. To counter this affect, 
the high frequency asymptote can be lowered by the same 
amount. Lowering the high frequency asymptote can be done 
by lowering the capacitive impedance, since, the capacitive 
impedance is inversely proportional to the capacitance, see the 
circuit in Figure B.l. Thus, this will also require multiplying 
the capacitance by the same proportional gain, K a . 

For the transverse approximation, illustrated in Figure B.3, 
it is seen that the high frequency asymptotes of the circuit 
approximation runs in parallel with the von Karman Spectral. 
Therefore, the asymptote corresponding to the transverse 



Figure B.2. — von Karman longitudinal spectral of Equation (B.2) 
and their respective circuit TF approximation of Equation (B.7), 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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Frequency, Hz 

Figure B.3. — von Karman transverse spectral of Equation (B.3) 
and their respective circuit TF approximation of Equation (B.7) 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Frequency, Hz 

Figure B.4. — Von Karman longitudinal acoustic wave spectral of 
Equation (B.2) and their circuit TF approximations with different 
magnitudes of adjustment (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


approximation needs to first be raised by doing the reverse of 
what was done to the capacitance in the previous step (i.e., 
multiplying the capacitance by the reciprocal of a certain 
correction factor, K c . Secondly, the magnitude of the 
transverse approximation can be adjusted upwards by a 
desired level, K a , by performing the two steps described for 
the adjustment of the longitudinal approximation. These 
correction factors can be computed as follows. 

K ta =10 AdB/2 ° 

(desired AdB level may range from 0 to 6 dB) 

Kl,c = 1, K VtC = S v y K (f c )-W v Jjc), 


(B. 15) 
(B.16) 


where f c is the frequency where this correction factor is to be 
computed. With these factors computed, Equation (B.13) for 
the capacitance can be modified utilizing Equations (B.15) and 
(B.16). But the approach would be to keep the values of the 
capacitance and resistance generic and instead incorporate 
these adjustment factors directly into Equation (B.7) as 


m,aA=- 


K, 

r K V 7 * 

^ t,a 

K t , c 


t,A 


-W t 


R t C t s 


+ 1 


(B.17) 


where 


K t9 A=K t9a Kf (B.18) 


and Kf (for the type of disturbance) comes from Equations 
(B.5) and (B.6) raised to the power r. 



Figure B.5. — Von Karman transverse acoustic wave spectral of 
Equation (B.3) and their circuit TF approximations with different 
magnitudes of adjustment (e=8.6e-5 m 2 /sec 3 , L=7 62 m). 


Figures B.4 and B.5 show the longitudinal and transverse 
acoustic wave spectral for the von Karman forms of Equations 
(B.2) and (B.3) and the equivalent circuit adjustments based on 
Equation (B.17) for different magnitude adjustments. The K u 
values for these adjustments are based on Equation (B.15). The 
Ky /C value used is 1.424. As discussed before, the important 
frequencies to match are the higher frequencies, above 1 Hz, 
because a typical control system design should have no problem 
attenuating disturbances at the lower frequencies. Therefore, the 
6 dB or even the 3 dB, adjustments would normally provide 
sufficient approximations to the von Karman spectral. 
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B.2 Temperature Disturbance 

Unlike the von Karman model forms of Equations (B.2) and 
(B.3) for the acoustic wave given by Soreide and Tank (1996), 
no comparable models exist for temperature, temperature 
generated acoustic wave, pressure, or density disturbances. 
Thus the approach here is to utilize the Kolmogorov spectral 
of Equation (B.l) for the acoustic wave disturbance and 
compute the differences from that and those of the 
Kolmogorov Spectral of Equation (B.l) for temperature. Then 
scale the respective von Karman models of Equations (B.2) 
and (B.3) by these differences to generate the von Karman 
models for temperature, pressure and density. 

Figure B.6 shows this scaling, with the resulting von 
Karman spectral for the temperature disturbance. This scaling 
can be done by graphically finding or by computing 
the difference in dB magnitude (at any frequency using 
Equation (B.l)) between the Kolmogorov temperature and the 
Kolmogorov longitudinal spectral, and then multiplying the 
von Karman longitudinal spectral by this scale factor 
(corresponding to the difference in dB magnitude) to come up 
with the von Karman temperature spectral. As a result, the von 
Karman temperature spectral is computed as 


Sr,VK{k) 

= 7.0s 2^5/3 1 , 

[l + (l.339(2Ti)z,q 2 ] 5/6 
((k 2 * m/sec )/H z) 


(B.l 9) 


The horizontal asymptote for this spectral, based on 
Equation (B.l 9), is 


K t - 14.0s 2/3 T 5/3 (B.20) 

The circuit approximation for this model is similar to that of 
Equation (B.l 7), except for some differences in the 
proportional gains and the fractional exponent as 


W TM =- 


K 


T,A 


-W T 


K T,a 
K T ,c 


RjCt^ 


+ 1 


(B.21) 


Where K Ta is the same as Equation (B.l 5) for the desired 
adjustment and 


Kt,a ~ Kt,ciKt > k t,c = 1 (B.22) 

As before, the fractional exponent is q=xr. However, this time 
for the temperature disturbance, the units conversion exponent 
is r=l/2 (see Fig. 1 for the temperature units). As such, the 


fractional exponent for temperature is q= 5/6, instead of q= 5/9 
for the acoustic wave disturbances. Figure B.7 shows a plot of 
the circuit approximation for temperature based on Equation 
(B.21), with K T>a = 1.4125 (i.e., raising the circuit 

approximation magnitude by 3 dB) , compared to its scaled 
von Karman spectral of Equation (B.l 9). 



Figure B.6. — Kolmogorov spectral of the longitudinal wave and 
temperature and scaling of the longitudinal von Karman 
spectral to come up with the von Karman temperature model 
form (e = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Figure B.7. — Kolmogorov and von Karman spectral for 
temperature and its von Karman circuit approximation 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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(B.25) 


By perturbing the relation of the speed of sound and 

temperature, a = yfjRT , and substituting the resulting 

expression for A a into the Mach number equation of M=v/a, 
the relation between the change in temperature and acoustic 
velocity can be obtained as 


My R Arri 

Av = — — AT (B.23) 

2 a 0 

Where y is the ratio of specific heats (-1.41), R is the universal 
gas constant (-287 - N*m/(kg*K)), and a Q is the standard 
speed of sound at that altitude (- 295 m/sec at supersonic 
cruise altitude). Figures B.8 and B.9 show the von Karman 
spectral for the longitudinal and transverse acoustic waves 
respectively, the temperature gust based on Equation (B.19) 
turned to either longitudinal or vertical acoustic wave velocity 
gust, and the combination of the two. As shown in these 
figures, for the combined acoustic wave the acoustic wave 
velocity due to temperature gust dominates in the low 
frequency range. At higher frequencies, both affect the total 
and the resultant slope is more of a combination of the two 
disturbances. The question in encountering atmospheric 
turbulence is if it’s possible that the two effects (i.e., the 
acoustic wave due to pure wind gusts and acoustic wave 
disturbance due to a temperature fluctuation) can combine at 
the same time to produce worst case conditions. As for 
longitudinal and transverse disturbances combining, in terms 
of the propulsion system, the assumption is that the 
longitudinal disturbance would be the worst case (i.e., the 
longitudinal and the transverse or vertical would not combine 
at the same time). This is especially true for propulsion 
systems located under the wing, since, the vertical acoustic 
wave loses some velocity when it’s turned to a longitudinal 
acoustic gust via the vehicle wing forebody (Ahsun (2004)). 


B.3 Pressure Disturbance 

For an atmospheric pressure disturbance, as discussed in the 
previous section, the Kolmogorov spectral of Equation (B.l) 
will be used to come up with the von Karman model form by 
scaling the respective models of the longitudinal acoustic 
wave disturbance. Figure B.10 shows this scaling, with the 
resulting von Karman spectral for the pressure disturbance. As 
a result, the von Karman pressure spectral is computed as 


Sp,vx(k) - 


5 8 s 2/3 Z 5/3 - 

[1 + (1.339(2ti)Mc) 2 ] 5/6 ’ 


(B.24) 


((Pa 2 *m/sec)/Hz) 

The horizontal asymptote for this spectral, based on 
Equation (B.24), is 


K P = 1 1.6s 2/3 Z 5/3 


Similar to the Temperature spectral, the circuit 
approximation for this model is the same as that of Equation 
(B.17), except for some differences in the proportional gains 
as 


W } 


K 


P,oA 


P,A 


-W P 


f v \ 1/ <1 
K P,a 

Kp, c 


RpCps 


+ 1 


(B.26) 



Figure B.8. — Von Karman longitudinal acoustic wave spectral, 
temperature disturbance turned to acoustic velocity, and the 
two combined. 



Figure B.9. — Von Karman transverse acoustic wave spectral, 
temperature disturbance turned to acoustic velocity, and the 
two combined. 
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Figure B.10. — Kolmogorov spectral of the longitudinal wave 
and pressure and scaling of the longitudinal von Karman 
spectral to come up with the von Karman pressure model 
form (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Frequency, Hz 

Figure B.1 1 . — Kolmogorov and von Karman spectral for 

pressure and its von Karman circuit approximation, (s = 8.6e-5 
m 2 /sec 3 , L = 762 m). 

Where K Pa is the same as Equation (B.15) for the desired 
adjustment and 


Kp,A - K P a K r p , K PtC = 1 (B.27) 

As with the temperature in the previous section, the fractional 
exponent is q=xr , with r=l/2. As such, the fractional exponent 
for pressure is the same as that for temperature (i.e., q= 5/6,). 
With the value of K Pa = 1.4125. Figure B.ll shows a plot 
of the circuit approximation for pressure based on 
Equation (B.26), compared to its scaled von Karman spectral 
of Equation (B.24). 



Frequency, Hz 

Figure B.1 2. — Kolmogorov and von Karman spectral for 
density and its von Karman circuit approximation (e = 8.6e-5 
m 2 /sec 3 , L = 762 m). 


B.4 Density Disturbance 

The density disturbance won’t be needed when acoustic 
velocity, temperature and pressure disturbances are included 
in a simulation. However for completeness, density 
disturbance will be discussed here in brief. The magnitude of 
the density disturbance can be computed from the equation of 
state, P=pRT, by perturbing this equation, which can be 
solved as 


P + a p p 

Ap = — +jl_ 

R(T a + AT) RT 0 


(B.28) 


where P Q and T 0 are the standard atmospheric pressure and 
temperature in Pa and Kelvin respectively, which would be 
around 5500 Pa and 216 K at supersonic cruise. Figure B.12 
shows a plot of the Kolmogorov, the Von Karman, and the 
circuit approximation spectral for density, all based on 
Equation (B.28) from values of temperature and pressure 
obtained in the previous sections. It can be seen from this 
figure that the Kolmogorov spectral obtained from Equation 
(B.28), actually looks like a finite von Karman spectral at low 
frequencies. Also, its low frequency asymptote would be 
approximately -20 dB, which equates to a density of 
0.1 kg/m 3 . From standard atmospheric tables, the atmospheric 
density at about 62,000 ft (i.e., at approximate cruise altitude) 
is also approximately 0.1 kg/m 3 . Even though, at the limit as 
frequency approaches zero, the Kolmogorov spectral 
approximates the mean atmospheric density at cruise, the 
appropriate density disturbance would be that shown for the 
von Karman or the circuit approximation. 
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Appendix C. — Fractional Order TF Fits of Atmospheric Disturbances 


Normally, with the derivation of circuit analogs and their 
respective TFs to describe atmospheric disturbances (as they 
were presented in the Appendix B), time domain simulations 
would be straightforward. However, these circuits and their 
TFs are fractional, which complicates performing time-domain 
simulations. The complication comes from a certain property 
of fractional order differential equations, which deals with 
non-locality. There are several works that deal with time 
domain solutions of fractional order differential equations 
Schmidt (2006) and Lorenzo (2008), but even to date, solving 
these types of differential equations is still challenging. Unlike 
integer order differential equations, the solution space of 
fractional order derivatives is not local and it depends on the 
whole prior history of the state, resulting in a state transition 
matrix that is in the form of a convolution integral. Therefore, 
solving the convolution at each time instant, taking into 
account all the prior history of the state, can be both 
complicated and taxing on the computing resources. 

In this section the attempt will be to formulate an integer 
order approximation of the fractional order TF. Given the 
parameters describing an atmospheric disturbance, like eddy 
dissipation rate and integral length scale, the formulations will 
explicitly solve for the integer order TF poles and zeros that 
approximate the fractional order disturbance. This approach 
avoids the laborious process of hand fitting such an 
approximation, for every type of disturbance, and for every 
time an atmospheric disturbance parameter changes. After 
forming this approximation to the fractional order TF, and by 
selecting the frequencies of interest for the disturbance as an 
input to this TF approximation, the time domain atmospheric 
disturbance can be formed. 

If the desire is to estimate the fractional order TF at each 
frequency decade by a single first order pole-zero pair, then 
the estimation process for each frequency decade starts with a 
pole followed by a zero (i.e., (o pi < (D zi ), with one last pole 
placed after the last decade in order to keep attenuating the TF 
magnitude as 


W Uo =K 


n 


s ! ®zi + 1 


. 1 s/(o vi +1 
V 1=1 / 


S ! ^pn + 1 1 


w t9 1 = 1,2,...,* (C.l) 



Figure C.l. — Pictorial diagram of fractional order TF fit. 


order TF midway through, in the first decade. A fractional 
order TF gain would be decreasing at a rate g*20 dB/decade, 
where q is the fractional order. Therefore, midway through the 
first decade the initial gain is expected to decrease by a factor 
of (V 2 ) *q *20 dB as 


K 

s/®pi + 1 


= ioT ?(20)/2 V 

.s-1 / 2decade 


(C.2) 


Divide by 20 is because of the 201og 10 (x) base scale. Solving 
for the frequency we get 



5=1/ 2 decade 


s 

lOd/2)^ _1 


(C.3) 


Also, midway in the first decade, the fractional TF, based on 
its circuit approximation in Appendix B, would have the same 
gain as that of Equation (C.2) as 


K 

(R,C t syi + \ 


s=l/ 2 decade 


-1 

= 10 ~2 q K 


(C.4) 


As was discussed in Appendix B, the natural frequency of the 
fractional order TF in Equation (C.4) is 


Where n is the number of pole-zero pairs of the estimated TF 
and K is the DC gain. So if the desire is to approximate the TF 
for three consecutive decades, with one pair of pole/zeros for 
each decade, starting at some desired frequency, then n would 
be three in this case. A sketch of this TF approximation as a 
staircase symmetrically located on top of the fractional order 
TF is shown in Figure C.l, with more description to follow 
later. 

Thus, if one first order pole-zero pair is used to estimate 
each decade of the fractional order TF, the first pole should be 
placed at a certain frequency such that the pole gain dropping 
at a rate of -20 dB/decade, intersects the gain of the fractional 


co 


n 


1 


(C.5) 


Knowing <n„, Equation (C.4) can be solved for s to find out at 
what frequency to place the first pole in order for the gains of 
these two TFs to intersect midway through the first decade as 


= CO„(l0( 1/2 )?-l) 1/<? 

s =1/2 decade 


(C.6) 


Substituting Equation (C.6) into Equation (C.3) for s, the 
frequency of the first pole can be solved as 
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CO^.1 = CO„(l0 (1/2) ^ — l) q 


(C.12) 


(C.7) 


a:(^/q } z1 +i) i 
s/<a pi +l l sa, “ tode 




So far in this development the assumption was that each 
decade of the fractional order TF will be approximated with 
one pole-zero pair as shown in Equation (C.l). Based on that, 
the estimated product of first order poles and zeros of 
Equation (C.l) will match the fractional order TF at half 
decade intervals. For that, the only distinction in 
Equation (C.7) is the 1/2 exponent in the right side of the 
equation. This successive ratio of desired decay matching, 
expressed here as H , can be generalized as 


<V =co„(l0( // /’>^ , -l)T (C.8) 

where H pX is the decay ratio for the first pole, where the two 
TFs would match or the magnitudes will intercept. 

Following the placement of the first pole, the need is to find 
out at what frequency to place the first zero of this pair. At the 
end of the first decade, the gain of the fractional TF would 
have dropped by 1 0“^ K as 

7 \s at first decade “ 1 ^ (C.9) 

(R t c t s) q +V 

By substituting 1/oy, using Equation (C.5), for R t C h the 
frequency s where the gain has dropped by the above 
magnitude can be calculated similarly as Equation (C.6) 

H = co„(l0^ -lj lq (C.10) 

at 1 decade 


Or in general terms, as a function of the second portion of the 
approximation, H zX desired to be matched (i.e., the first 
portion at Vi decade was matched with a pole, the second 
portion at 1 decade is matched with a zero (H zX = 1), and so on), 
this frequency can be expressed as 

co/fei=®„(l0(^.)9-l J /q CC.11) 


where 


From this equation, the frequency of the first zero can be 
calculated as 


©zl = 


©ifel 

\Q~( H :') q (m Hz \/m p \ +l)-l 


(C.13) 


Thus, choosing the ratio of decades at which it’s desirable to 
match the TF estimate, Equations (C.2), (C.8), (C.ll), and 
(C.13) can be used to compute the frequency of the first zero. 

Following the same procedure to calculate the frequency of 
the second pole, by using Equation (C.9), but with the gain 
lO~ (V2)q K, (i.e., H p 2=1.5 decades), 

(C.14) 


Then using the TF with two poles and one zero equated to this 
gain, similar to Equation (C.12), the frequency of the second 
pole can be calculated as 


©/?2 — 


^Hpli^Hpl/^pl +l) 

lO^Mo^/cOzi + l)- 1 


(C.l 5) 


The same procedure can be followed to calculate the rest of 
the frequencies of the poles and zeros of the estimating TF. In 
order to generalize these formulations (for approximating the 
fractional order TF), let’s define p pz as the desired density of 
pole-zero pairs per frequency decade. In addition, let’s also 
define the number of equal decade subintervals (in log 
frequency scale) as r|, where a pole or a zero will be used to 
closely match this TF as 


1 

2p P z 


(C.l 6) 


Then, pole and zero vectors of decade intervals can be formed, 
where these TF will be closely matched as 


(decade) 

The purpose of the co H frequencies is to ensure symmetry (i.e., 
to symmetrically locate the TF approximation on top of the 
fractional order TF. Thus, in order to accomplish symmetry, 
the desire is to also have the estimating TF with the first pole- 
zero pair to equate to the same gain at this frequency, co H zi, 
(i.e., the same gain as the right hand side of Equation (C.9), 
see also Figure C.l, as 


H p = j\ [2(1) - 1 2(2) - 1...2 (m p ) - 1] (C.17) 

H z = ^[2(l)2(2)...2(m z )\ (C.18) 

For n number of decades over which to estimate the fractional 
order TF, m p and m z can be computed as 

m p = (n- l)/r| (C.l 9) 

m z = m p - 1 (C.20) 
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For instance, for a density of one pair of pole-zeros per decade 
(i.e., p pz =1) and estimating the TF for three decades (i.e., 
n=3), these vectors would have the values H pi = [1/2 3/2 5/2 
7/2], H zi = [1 2 3]. 

Inspecting Equation (C.17), it can be seen that the first 
element of the vector, corresponding to the calculation of the 
first pole frequency, will always be equal to the value of r\. 
Therefore, for convenience, Equation (C.8) can also be 
expressed in terms of r|. 


<D / , 1 =(D„(lO’W-l) 1 / (C.21) 

Also, inspection of the equations that describe the 
frequencies of the zeros and poles (i.e., Eqs. (C.ll), (C.13), 
(C.14) and (C.15) respectively, it can be seen that a 
generalized relationship for these equations can be written as 
follows 


0:1 Hpi = ®« (l -\ ) (C.22) 

A more generalized expression of Equation (C.22) can be 
arrived by substituting into it H p from Equation (C.17) in 
terms of the pole number, i and rj as 

(0 Hpi = ©«(10 n(2M)? -l) l/</ i = 2 (C.23) 

Then a generalized expression for the calculation of the pole 
placement frequencies can be formulated as 

= 

_ Kg pi® Hpi {®Hpi / 03 pi - 1 + 1 Jto H pj / pj-2 + 1) ■ • • {to Hpi / 03 pil + -Q 
IQ ( H pi) q ]^CO Hpi j CO zi _\ + l)(c 0 Hpi / (D zi _ 2 + l)...(cO// p i- / (D z i + l)]- 1 

\ (C.24) 

K apitoHpi± \ \^Hpi /topi-j +1) 

j=\ . 0 

= tzt > i = 2,...,m p 

io«n(^M ; , y+ i)-i 

7=1 

The proportionality factors K api have been inserted in case any 
final adjustments are needed to these frequencies to improve 
the TF fit. 

Similarly, for the frequencies for placing the zeros 

(o Hzi =(o„{loM-lf q i = l,2,...,m z (C.25) 

Or by substituting in Equation (C.21) the H z vector from 
Equation (C.18) in terms of the zero number, z, and r\ as 

®^=®„(l0 2 w-l f q i = 1,2, ...,m z (C.26) 

Then a generalized expression for the calculation of the zero 
placement frequencies can be formulated as 


tozi = 

1 0 ( H zi)q [{a Hzi /(D pi + l)((0 Hzi l<Q pi -\ + l) .. i^Hzi /topi + 1)] ■ - 1 
i - 1 

K(ozi ® Hzi Hzi f^zi-j "l" 1 ) 

= — for / = l,2,...m z 

i 

10- 2 w^Q( C0ffz ./ tO;) . + i)_i 

i = i 


(C.27) 


These formulations complete the computation of the pole-zero 
frequencies for the TF approximation of fractional order 
atmospheric disturbances. In terms of the indexes defined 
here, Equation (C.l) can be rewritten as 

m z 

| ( s/®zi l) 

W t , 0 = K tJlt -J W t (C-28) 

J m p 

n(v®/rf + i) 

i 

As described in Appendix B, the von Karman spectral, 
Equations (B.2) and (B.3), can be approximated with a 
fractional order circuit with a natural frequency that can be 
determined using Equation (B.8). Equation (B.8) is 
conveniently reproduced below in Equation (C.29). 

o\=-^- (c.29) 


with an equivalent capacitance and resistance for the type of 
disturbance t , derived as Equation (B.13) and (B.14) 

C= I 

(a t z 2l3 J lx (2nMa) 

R t = \339(2n)(a t & 2 / 3 f x L 

Based on these derivations, TF approximations of fractional 
order atmospheric turbulences can be developed (for the von 
Karman spectral, for any fractional order o < <7 < 1) by using 
Equation (C.28), together with supporting Equations (C.l 6), 
(C.l 9) to (C.21), (C.25) to (C.27), (C.23) to (C.24) and (C.26) 
to (C.27), (C.29) to (C.31). 

In a pictorial sense, by referring back to Figure C.l, the first 
order pole-zero pair TF approximation of a fractional order 
atmospheric turbulence can be viewed as a staircase 
symmetrically located on top of the fractional order TF. The 
utility of (o Hpi and co Hzi are to maintain this symmetry. This is 
accomplished in the derivations above by formulating the 
symmetry frequencies, CD/fs, such that the length segments 
between intercepts (like d 1 and c/2) or alternatively the 
horizontal distances on either side of the intercepts are equal. 

In the subsections that follow, more detail will be discussed 
about developing the approximations for each type of 


(C.30) 

(C.31) 
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disturbance. In addition, some final adjustments that need to 
be made to the pole-zero frequencies and the proportionality 
constants K will be presented. Without these final adjustments, 
the fits will resemble the example fit shown in Figure C.2, 
with the fractional order fit winding its way around the 
fractional order TF 


C.l Longitudinal and Transverse Acoustic 
Wave Disturbances — TF Fits 


By applying the formulations derived in Section C.l, TF fits 
for the longitudinal and transverse atmospheric disturbances 
can be constructed based on Equation (C.28). Figure C.3 
shows a plot of the longitudinal von Karman disturbance of 
Equation (B.2), it’s circuit TF approximation reproduced from 
Figure B.2, and the TF fit computed using the equations 
derived previously in this Appendix for q=xr= 5/9 (i.e., x=5/3, 
7=1/3), r| = l/2, n= 3 (decades), and with atmospheric 
turbulence parameters s=8.6e-5 (m 2 /sec 3 ) andZ=762 m. 

The maximum error in this TF fit for the circuit 
approximation, using 4 poles and 3 zeros, over a span of 3 l A 
decades, amounts to approximately 1.5 dB. The error for the 
actual von Karman spectral is larger; about 4 dB at 0.1 Hz and 
decreases as the frequency increases. Assuming, that the 
control system design can sufficiently attenuate frequencies at 
this low frequency range (as it should) then this TF fit 
approximation can be acceptable for the actual von Karman 
spectral. For the TF fit shown in Figure C.3, which is 
an approximation to the longitudinal von Karman spectral 
of Equation (B.2) the pole and zero frequencies in 
Equation (C.28) are 


(Opt = [0.6 12.54 85.71 2062.0], 

r .. (rad/sec) 

©*■=[3.82 22.92 312.38] 


(C.32) 


with K /flt in Equation (C.28) given by the numerator in 
Equation (5) raised to 1/3 power to account for the units 
conversion as 


K, m =(5.4fi 2/ V 5/3 )3 (C.33) 

The preceding development to derive the equations to 
approximate a fractional order TF was based on a single 
fractional order TF, and a circuit approximation to the 
atmospheric disturbance was used for this development. 
However, as indicated by Equation (B.2) and (B.3), the von 
Karman spectral is not that of a single order fractional TF, 
especially in the lower frequency spectral. Therefore, for 
better accuracy, the derived equations can be adjusted to better 
fit the von Karman spectral. Examining Figure C.3, it can be 
seen that the TF fit can better approximate the von Karman 



Frequency, Hz 

Figure C.2. — An example Fit without final adjustments to 
proportionality factors and frequencies. 



Frequency, Hz 

Figure C.3. — Longitudinal von Karman spectral, circuit 
approximation, and TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


spectral if the natural frequency co„, Equation (C.5), is 
increased, about proportionally to the difference in magnitudes 
between the fractional order TF and the approximation TF, at 
that frequency. However, doing that will also increase the 
pole-zero frequencies of the TF fit based on the equations 
derived. Therefore, some equivalent adjustment to some of 
these frequencies will need to be made. Figures C.4 and C.5 
show the TF fit for the longitudinal and transverse von 
Karman disturbances respectively, with the following 
adjustments to Equations (C.24), (C.29) and (C.27) 
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K-l / v,co — \^am 9 K($pi , K^ z i ] 

= [2.4; 1 1 1/2.4 1/1.5; 1 1 l] 


(C.34) 


With Kifaco symbolizing the proportionality constants for the 
frequency adjustments for both the longitudinal and transverse 
disturbances, which gives the frequency values of ( for e=8.6e- 
5 m 2 /sec 3 and L = 762 m) 


(o p i = [l .46 30.10 85.71 1593. l], 
oo z/ = [9.18 55.02 335.48] 


(C.35) 


The proportionality factor, Kj >rit , in Equation (C.28) is given 
again by Equation (C.33). Similarly, K v> f lt is developed using 
Equation (B.6) as follows, but also adjusting by 3 dB which 



Figure C.4. — Longitudinal von Karman spectral, adjusted TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


will appropriately raise its magnitude in order to match the 
high frequency asymptote as 

K v ,fit = 1.4(2. 7e 2/3 Z 5/3 )j (C.36) 

As can be seen in Figure C.5, the TF fit of the transverse 
disturbance is not very accurate at the low frequency, which 
would normally be fine for control system design purposes as 
discussed before. But it may not be acceptable, for let’s say, 
wing loading. Thus, further adjustments could be made if 
desirable. For instance, instead of multiplying K v>& t with the 
3 dB gain of 1 .4, the natural frequency cd„ of Equation (C.29) 
could be multiplied by this factor to raise its frequency. Then 
for final adjustments, the natural frequency can be also 
multiplied by the ratio of the difference of frequencies between 
the fit and the actual von Karman spectral obtained at some 
fixed amplitude, at some high frequency. This can be done 
either graphically or analytically. Performing these 
manipulations, a better fit of the transverse disturbance is shown 
in Figure C.6, with a final value of K om = 4.27. This will result 
in the following adjustments for the transverse disturbance 

*v,a> — [^00/7 ? ^Gipi ? K 0 , zl ] ^ ^ 

= [4.27; 1 1 1/2.4 1/1.5; 1 1 l] 

which gives the frequency values of ( for e=8.6e-5 m 2 /sec 3 and 
L=162 m) 


(0 ViP i = [2. 60 53.56 152.55 2835.3], 
(0 VZ i = [l6.33 97.92 597.07] 

with a K v f [t expression as 

K v ,m={2Je 2/3 L 5/3 )\ 


(C.38) 


(C.39) 




Figure C.5. — Transverse von Karman spectral, adjusted TF fit Figure C.6. — Transverse von Karman spectral, final adjusted 

(s = 8.6e-5 m 2 /sec 3 , L = 762 m). TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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C.2 Temperature Disturbance — TF Fit 

As discussed before, the units conversion factor for 
temperature is r= 1 / 2 , which makes the fractional exponent 
value g=5/6 (x=5/3). Substituting this value for the fractional 
exponent, with the same values of r|, n, and with the same 
atmospheric parameter values as those used for the 
longitudinal and transverse acoustic disturbances, and again 
using the formulations previously developed in this Appendix, 
the poles and zeros for the TF fit can be directly computed as 
before. The resulting TF fit with K t j x{ in Equation (C.28) given 
by Equation (B.20) raised to the l A power as 

K t =( 14.0s 2 7 V 5/ A (C.40) 

This fit (not shown here) will end-up with the same low and 
high frequency asymptotes as with the von Karman form, but 
at somewhat lower frequency. Similarly to the longitudinal 
disturbance, this can be corrected by multiplying K ro „ in 
Equation (C.29) by 3/2, which turns out to be the ratio of the 
frequencies of these two spectra at some high frequency, for 
select fixed amplitude. Finally, by some minor adjustments to 
the high frequency poles, these spectra can be matched fairly 
close, as shown in Figure C.7. This will result in the following 
adjustments for the temperature disturbance 

K-T, go — Kn ? Kcopi j Kazi ] (C 4 1 ) 

= [1.5; 1 1 1/1.1 1/1.2; 1 1 l] 

which gives the frequency values of (for s=8.6e-5 m 2 /sec 3 and 
L= 762 m) 


co 7 1 D i = [l . 1 0 25.11 109.77 816.351 

(C 42) 

©r,zi= [33.04 45.64 602.36] 
with a K t flt , using Equation (B.20) as 

*7\fit(temp) = (l4.0s 2 /3L5/3)1 /2 (C.43) 

If the desire is to simulate the acoustic disturbance 
generated by a temperature gust, this conversion can be done 
by utilizing the perturbation relation of the speed of sound 
with temperature described in Appendix B (Eq. (B.23)) as 


Av = 


MyR 

2a 0 


AT 


(C.44) 


Together with Equation (C.36) to adjust K T f lt as 

£r,fit(acoustic) =^Vl4.0s 2/3 P /3 (C.45) 

zu 0 



Figure C.7. — Temperature von Karman spectral and its TF fit 
(s=8.6e-5 m 2 /sec 3 , L=7 62 m). 



Frequency, Hz 

Figure C.8. — Von Karman acoustic wave velocity spectral due to 
temperature gust and its TF fit (c = 8.6e-5 m 2 /sec 3 , L = 762 m). 


The acoustic velocity disturbance due to a temperature gust, at 
approximate supersonic cruise altitude is shown in Figure C. 8 . 

Based on all the derivations carried out so far, the TF fits 
should work for different values of atmospheric parameters 
(like those of Eq. (C.45)). Figures C.9 and C.10 as compared 
to Figures 3 and 4 for their shape, show the TF fits with 
different values of eddy dissipation rates and integral scale 
lengths, by employing the same TF fit equations and without 
changing any other parameters, like the K T >G) values, 
Equation (C.41). 
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to temperature gust and its TF fits for different integral scale 
lengths (s = 8.6e-5 m 2 /sec 3 ). 



Figure C.1 1 . — Pressure von Karman spectral and its TF fit 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Frequency, Hz 

Figure C.10. — Von Karman acoustic wave velocity spectral due 
to temperature gust and its TF fits for different eddy 
dissipation rates (L=7 62 m). 


(C.27), (C.23) to (C.24) and (C.26) to (C.27), (C.29) to (C.31), 
the poles and zeros for the TF fit for pressure can be computed 
as before. The resulting TF fit, with K t f lt in Equation (C.28) 
given by Equation (B.25) raised to the l A power, is shown in 
Figure C.ll. It turns out, that adjustments to K P >Q) are exactly 
the same as those for K T}(0 , Equation (C.41). For completeness, 
these frequency adjustments and the K P f lt equation for 
pressure disturbances are listed below. 

Kp 9 c o = i K($pi , K^ z f ] 

= [l.5; 1 1 1/1.1 1/1.2; 1 1 l] 

which gives the frequency values of (for s=8.6e-5 m 2 /sec 3 and 
L=162 m) 

(Op pi = [l.lO 25.11 109.77 816.35], 

(Qp : zi = [33.04 45.64 602.36] ^ ’ 

with a Kp'fa , using Equation (B.25) as 

K P ' flt = Vll.6s 2/3 L 5/3 (C.48) 


C.3 Pressure Disturbance — TF Fit 

The units conversion factor for pressure is the same as that 
for temperature (i.e., r=l/2), which also makes the fractional 
exponent value q=5/6. Substituting this value for the fractional 
exponent, with the same values of r|, n, and H as those used 
for the longitudinal, transverse and temperature disturbances, 
and again using Equations (C.16), (C.19) to (C.21), (C.25) to 


C.4 Density Disturbance — TF Fit 

As discussed in Section B.4, once the temperature and 
pressure TF fits have been calculated based on the 
formulations in the previous two sections, the density 
disturbance can be directly obtained by utilizing 
Equation (B.28). A plot of this density spectral is shown in 
Figure C.12. 
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Figure C.12. — Density von Karman spectral and its TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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